﻿#include "discrete_face_torus.hxx"
#include "discrete_face_loop.hxx"
#include "idx_mesh_utils.hxx"

/**
 * @brief 初始化环面对象及其相关参数
 *
 */
void torus_faceter::init() {
    // 将f对象中的几何形状更新方程的引用转换为环面类型的指针
    // setupTorus(f);
    tori = (torus*)&f->geometry()->equation_for_update();
    // 获取参数u和v的取值范围
    SPAinterval u_interval = f->geometry()->equation().param_range_u();
    SPAinterval v_interval = f->geometry()->equation().param_range_v();
    // 参数u的起始点和结束点
    u_le = u_interval.start_pt();
    u_ri = u_interval.end_pt();
    // 参数v的起始点和结束点
    v_le = v_interval.start_pt();
    v_ri = v_interval.end_pt();
    // 设置容差值，用于判断是否跨越边界
    u_tolerance = SPAresabs;
    v_tolerance = SPAresabs;
    // 如果在u(v)方向上是封闭的，将容差值设为u(v)参数范围的一半
    if (f->geometry()->equation().closed_u()) u_tolerance = tori->param_range_u().length() / 2;
    if (f->geometry()->equation().closed_v()) v_tolerance = tori->param_range_v().length() / 2;

    // 设置左下角点的坐标
    le_low.x = u_le;
    le_low.y = v_le;
    // 设置右上角点的坐标
    ri_high.x = u_ri;
    ri_high.y = v_ri;

    SPAunit_vector zeroNormal = f->geometry()->equation().eval_normal(SPApar_pos(0, 0));
    if (f->sense()) zeroNormal = -zeroNormal;
    if (zeroNormal != tori->uv_oridir) normalR = REVERSED;
}

/**
 * @brief 计算u和v方向上的步长
 *
 */
void torus_faceter::decideUVlen() {
    // 返回ulen和vlen表示u和v方向上的分割数量
    // get_uvlen_torus(f, ulen, vlen, nt, st);
    decideUV(f, ulen, vlen, nt, st);
    // 计算u和v方向上的步长
    du = (u_ri - u_le) / ulen;
    dv = (v_ri - v_le) / vlen;
    // printf("%lf %lf\n", du, dv);
    // printf("====================%d %d==================\n", ulen, vlen);
}

// 满足表面容差的公式
/*
    cos(b/2)(R+|r|cos(a/2))≥R+|r|-st
    令t=R+|r|-st，其中a,b分别为u和v方向上的步进值。
    令a=2b,代入上述不等式，解方程得result=cos(b/2)，即b=2*acos(result)。
 */
double f_st(double R, double r, double t) {
    double inside_sqrt = (2 * std::pow(R, 3) - 6 * std::pow(R, 2) * r + 6 * R * std::pow(r, 2) - 2 * std::pow(r, 3) + 27 * r * std::pow(t, 2)) / r;
    assert(inside_sqrt >= 0);
    double sqrt_part = 6 * std::sqrt(3) * std::sqrt(inside_sqrt) + 54 * t;
    assert(sqrt_part != 0);
    double cube_root_part = std::cbrt(sqrt_part * std::pow(r, 2));
    assert(cube_root_part != 0);
    double result = cube_root_part / (6 * r) - (R - r) / cube_root_part;

    return result;
}

// 满足角度容差的公式
/*
    1.cos^2(a/2)(cos(b）-1)≥cos(nt*π/180.0）-1；
    2.cos^2(a/2)(cos(b）+1)≥cos(nt*π/180.0）+1； 当该不等式满足时，1号不等式一定满足。
    3.a≤nt*π/180.0
    令t=cos(nt*π/180.0）+1，其中a,b分别为u和v方向上的步进值。
    令a=2b,代入上述2号不等式，解方程得result=cos(b),即b=acos(result)。
*/
double f_nt(double t) {
    double part1 = -8 + 108 * t + 12 * std::sqrt(81 * t * t - 12 * t);
    double root1 = std::cbrt(part1);
    assert(root1 != 0);
    double result = root1 / 6 + 2 / (3 * root1) - 1.0 / 3;
    return result;
}

void torus_faceter::decideUV(FACE* f, int& ulen, int& vlen, double nt, double st) {
    torus* tori = (torus*)&f->geometry()->equation();

    double r = tori->minor_radius;
    double abs_r = std::abs(tori->minor_radius);
    double R = tori->major_radius;
    double abs_R = std::abs(tori->major_radius);

    double val = f_st(R, abs_r, R + abs_r - st);
    val = std::min(1.0, std::max(-1.0, val));
    double step_v1 = 2 * acos(val);
    double step_u1 = 2 * step_v1;

    double temp = f_nt(1 + cos(nt * M_PI / 180.0));
    temp = std::min(1.0, std::max(-1.0, temp));
    double step_v2 = acos(temp);
    double step_u2 = 2 * step_v2;

    double step_utemp = std::min(step_u1, step_u2);
    double step_u = std::min(step_utemp, nt * M_PI / 180.0);  // 需满足角度容差的3号不等式
    double step_v = std::min(step_v1, step_v2);

    ulen = static_cast<int>(std::ceil(2 * M_PI / step_u));
    vlen = static_cast<int>(std::ceil(2 * M_PI / step_v));

    // @todo 确定最小分段数
    int min_ulen = 3;
    int min_vlen = 3;
    ulen = std::max(ulen, min_ulen);
    vlen = std::max(vlen, min_vlen);
}

bool torus_faceter::getLoops() {
    return false;  // 返回false
}

void torus_faceter::peripheryProcess() {
}

bool torus_faceter::USeperationProcess() {
    return false;  // 返回false
}

bool torus_faceter::VSeperationProcess() {
    return false;  // 返回false
}

void torus_faceter::holeProcess() {
}

/**
 * 存储扫描线与边界交点的信息
 */
struct scan_inter {
    double u;       // 扫描线与边界交点的 u 坐标值
    int low_inter;  // 扫描线与边界的交点中位于扫描线下方的交点数量
    int up_inter;   // 扫描线与边界的交点中位于扫描线上方的交点数量
};

void torus_faceter::performTriangulationAndAttach() {
    REVBIT r = f->sense();
    // 调用CDT
    CDT::Triangulation<double>* pcdt = nullptr;  // CDT进行三角剖分
    int _count = 0;
    // do {
    if (pcdt) ACIS_DELETE pcdt;
    pcdt = ACIS_NEW CDT::Triangulation<double>;         // 新建一个三角剖分对象。
    CDT::RemoveDuplicatesAndRemapEdges(points, edges);  // 处理点和边的重复，重新映射边的索引

    CDT::EdgeUSet markedEdge;
    for (int i = 0; i < loop_edge_siz; i++) {
        markedEdge.insert({ edges[i].v1() + 3, edges[i].v2() + 3 });
    }
    // 将处理过的点集和边集插入pcdt对象中
    pcdt->insertVertices(points);
    pcdt->insertEdges(edges);

    auto&& mp = CDT::EdgeToPiecesMapping(pcdt->pieceToOriginals);
    for (auto& x : mp) {
        if (markedEdge.find(x.first) != markedEdge.end()) {
            for (auto& v : x.second) {
                markedEdge.insert(v);
            }
        }
    }

    pcdt->eraseOuterTrianglesAndHoles(markedEdge);           // 清除外部三角形和孔洞
    CDT::TriangleVec triangles = pcdt->triangles;            // 获取三角形的集合
    std::vector<CDT::V2d<double>> vertice = pcdt->vertices;  // 获取顶点的集合
    std::vector<UV_VERTEX> uv_vertices;
    // 遍历每个顶点
    for (int i = 0; i < vertice.size(); i++) {
        CDT::V2d<double> pp = vertice[i];
        SPApar_pos param(pp.x, pp.y);  // 获取参数位置
        SPAposition pos = f->geometry()->equation().eval_position(param);
        SPAunit_vector normal = f->geometry()->equation().eval_normal(param);
        if (tori->degenerate() && fabs(param.u - u_ri) < SPAresabs) {
            param.u -= 2 * SPAresabs;
            normal = f->geometry()->equation().eval_normal(param);
            param.u += 2 * SPAresabs;
        }
        else if (tori->degenerate() && fabs(param.u - u_le) < SPAresabs) {
            param.u += 2 * SPAresabs;
            normal = f->geometry()->equation().eval_normal(param);
            param.u -= 2 * SPAresabs;
        }
        if (r) normal = -normal;
        // 将顶点信息添加到mesh对象中
        uv_vertices.emplace_back(pos, normal, param);
    }
    NewIndexedMesh(f, uv_vertices, triangles);
}

void removeOddEvenIntersections(std::vector<scan_inter>& inter_u) {
    for (unsigned i = 0; i < inter_u.size(); i++) {
        bool isOdd = (inter_u[i].up_inter + inter_u[i].low_inter) % 2 == 1;
        bool isEven = inter_u[i].low_inter % 2 == 0 && inter_u[i].up_inter % 2 == 0;

        if (isOdd && (i == 0 || i == inter_u.size() - 1)) {
            inter_u.erase(inter_u.begin() + i);
            i--;
        }
        else if (isEven) {
            inter_u.erase(inter_u.begin() + i);
            i--;
        }
    }
}

/**
 * @brief 根据给定的参数和边界条件在环面上进行扫描和交点计算，然后使用CDT进行三角剖分，生成索引化的网格，并附加到指定的面对象上
 *
 */
void torus_faceter::attachMesh() {
    SPApar_pos param;
    loop_edge_siz = edges.size();
    for (double v = v_le; v < v_ri + SPAresabs / 2; v += dv) {
        // 构造一条从y = v 的直线 求交
        // 求交结果记录在inters中，以两个点为一组，两个点中间的点为inside的点，添加到points中
        SPAposition a(u_le - 2 * SPAresabs, v, 0.0);
        SPAunit_vector dir(1.0, 0.0, 0.0);                         // x轴正方向
        straight v_line(a, dir);                                   // 从点a出发、方向为dir的直线
        SPAinterval line_range(0.0, u_ri - u_le + 4 * SPAresabs);  // 创建直线v_line的范围
        v_line.limit(line_range);
        // 用于存储扫描线与曲面的交点信息
        std::vector<scan_inter> inter_u;

        for (unsigned i = 0; i < edges.size(); i++) {  // 遍历edges向量中的每条边
            // 获取当前边的起始点和终点
            point curr = points[edges[i].v1()];
            point next = points[edges[i].v2()];
            // 跳过长度为零的边
            if (curr == next) continue;

            SPAposition b(curr.x, curr.y, 0.0);
            SPAposition c(next.x, next.y, 0.0);

            SPAvector edge_dir = c - b;  // 从b到c的方向向量
            straight edge_line(b, normalise(edge_dir));
            SPAinterval edge_range(0.0, (edge_dir).len());
            edge_line.limit(edge_range);
            curve_curve_int* results_head = int_cur_cur(v_line, edge_line);
            curve_curve_int* results = results_head;

            // 检查是否有曲线与曲线的交点
            if (results != nullptr) {
                // 如果有多个交点，就删除曲线与曲线的交点链表，继续下一个边的处理
                if (results->next != nullptr) {
                    delete_curve_curve_ints(results_head);
                    continue;
                }
                double p1 = results->param1;
                SPAposition inters = v_line.eval_position(p1);

                unsigned j = 0;
                for (; j < inter_u.size(); j++) {
                    if (inter_u[j].u == inters.x()) {
                        break;
                    }
                }
                if (j == inter_u.size()) {
                    // 将一个新的scan_inter结构体对象添加到inter_u向量中
                    inter_u.push_back(scan_inter(inters.x(), 0, 0));
                }
                // 如果交点在边界上方或下方
                if (is_equal(v, curr.y) && is_equal(v, next.y)) {
                    inter_u[j].low_inter += 1;                                                 // 位于下方的交点数量增加
                    inter_u[j].up_inter += 1;                                                  // 位于上方的交点数量增加
                }
                else if (is_equal(v, curr.y) && v > next.y || v > curr.y && is_equal(v, next.y)) {  // 如果交点在边界上方
                    inter_u[j].low_inter += 1;
                }
                else if (is_equal(v, curr.y) && v < next.y || v < curr.y && is_equal(v, next.y)) {  // 如果交点在边界下方
                    inter_u[j].up_inter += 1;
                }
                else {  // 交点同时在边界的上方和下方
                    inter_u[j].low_inter += 1;
                    inter_u[j].up_inter += 1;
                }
            }
            // 释放results_head指向的曲线与曲线的交点链表内存

            delete_curve_curve_ints(results_head);
        }

        // 对结构体对象按照u值进行升序排序
        std::sort(inter_u.begin(), inter_u.end(), [](const scan_inter& a, const scan_inter& b) { return a.u < b.u; });
        // 删除奇数次交点或偶数次交点
        removeOddEvenIntersections(inter_u);
        // 在指定的扫描线上生成新的点坐标，并将这些点坐标添加到points向量中
        for (int i = 0; i < inter_u.size(); i += 2) {
            double inf = inter_u[i].u;  // 获取当前段的起始u坐标
            if (i + 1 >= inter_u.size()) continue;
            double sup = inter_u[i + 1].u;  // 获取当前段的结束u坐标
            bool first = true;
            for (double u = get_cloest_boundary_point(u_le, du, inf); u < sup; u += du) {
                if (u < inf) continue;
                points.push_back(point(u, v));
                if (!first) edges.emplace_back(points.size() - 1, points.size() - 2);
                if (first) first = false;
            }
        }
    }

    // 进行三角剖分并附加到面对象
    performTriangulationAndAttach();
}

std::vector<SPApar_pos> generatePoints(SPApar_pos begin, SPApar_pos end, double v_le, double u_le, double dv, double du) {
    double min_x = std::min(begin.u, end.u);
    double max_x = std::max(begin.u, end.u);
    double min_y = std::min(begin.v, end.v);
    double max_y = std::max(begin.v, end.v);
    bool is_u = begin.u == end.u;
    std::vector<SPApar_pos> tmp;
    SPApar_pos tmp_point;

    if (is_u) {
        int k = static_cast<int>(std::floor((min_y - v_le) / dv));
        tmp_point.u = min_x;
        tmp_point.v = v_le + k * dv;
        for (; tmp_point.v <= max_y; tmp_point.v += dv) {
            if (tmp_point.u >= min_x && tmp_point.u <= max_x && tmp_point.v > min_y && tmp_point.v < max_y) {
                tmp.emplace_back(tmp_point.u, tmp_point.v);
            }
        }
    }
    else {
        int k = static_cast<int>(std::floor((min_x - u_le) / du));
        tmp_point.v = min_y;
        tmp_point.u = u_le + k * du;
        for (; tmp_point.u <= max_x; tmp_point.u += du) {
            if (tmp_point.u > min_x && tmp_point.u < max_x && tmp_point.v >= min_y && tmp_point.v <= max_y) {
                tmp.emplace_back(tmp_point.u, tmp_point.v);
            }
        }
    }

    // 根据排序条件排序
    bool is_less = is_u ? (begin.v < end.v) : (begin.u < end.u);
    if (is_u) {
        std::sort(tmp.begin(), tmp.end(), [&](const SPApar_pos& a, const SPApar_pos& b) { return is_less ? a.v < b.v : a.v > b.v; });
    }
    else {
        std::sort(tmp.begin(), tmp.end(), [&](const SPApar_pos& a, const SPApar_pos& b) { return is_less ? a.u < b.u : a.u > b.u; });
    }

    return tmp;
}

/**
 * @brief 处理不同的边界和环，对环面进行网格化处理，并将生成的网格附加到面对象上
 *
 * @return 处理结果
 */
outcome torus_faceter::facet() {
    init();  // 初始化
    decideUVlen();

    auto generator = [&](SPApar_pos begin, SPApar_pos end) { return generatePoints(begin, end, v_le, u_le, dv, du); };

    auto singularitygenerator = [&](SPApar_pos begin, SPApar_pos end) { return generatePoints(begin, end, v_le, u_le, dv, du); };

    auto loop_processor = LoopProcessor(f);
    loop_processor.RegisterGenerator(generator);
    loop_processor.RegisterSingularityGenerator(singularitygenerator);
    loop_processor.Execute();
    for (auto x : loop_processor.points) {
        points.emplace_back(x.u, x.v);
    }
    for (int loop_idx = 0; loop_idx < loop_processor.loops_indexs.size(); loop_idx++) {
        auto& point_indexs = loop_processor.loops_indexs[loop_idx];
        int point_indexs_size = point_indexs.size();
        if (point_indexs_size < 2) continue;
        for (int i = 1; i < point_indexs_size; i++) {
            edges.emplace_back(point_indexs[i], point_indexs[i - 1]);
        }
        edges.emplace_back(point_indexs.front(), point_indexs.back());
    }

    // 调用attachMesh()函数生成网格并附加到面对象上
    attachMesh();
    return outcome();
}

/**
 * @brief 对环面进行离散
 *
 * @param f 待处理的环面
 * @param nt 法向容差
 * @param st 距离容差
 * @return 返回网格化处理结果
 */
outcome facet_face_torus(FACE* f, double nt, double st) {
    torus_faceter faceter(f, nt, st);
    return faceter.facet();
}